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Two  approaches  are  commonly  used  for  handling  fric¬ 
tional  contact  within  the  framework  of  the  Discrete  Ele¬ 
ment  Method  (DEM).  One  relies  on  the  Complementarity 
Method  (CM)  to  enforce  a  non-penetration  condition  and 
the  Coulomb  dry-friction  model  at  the  interface  between  two 
bodies  in  mutual  contact.  The  second  approach,  called  the 
Penalty  Method  (PM),  invokes  an  elasticity  argument  to  pro¬ 
duce  a  frictional  contact  force  that  factors  in  the  local  de¬ 
formation  and  relative  motion  of  the  bodies  in  contact.  We 
give  a  brief  presentation  of  a  DEM-PM  contact  model  that 
includes  multi- time -step  tangential  contact  displacement  his¬ 
tory.  We  show  that  its  implementation  in  an  open  source  sim¬ 
ulation  capability  called  Chrono  is  capable  of  accurately 
reproducing  results  from  physical  tests  typical  of  the  field  of 
geomechanics;  i.e.,  a  direct  shear  tests  on  a  mono-disperse 
material.  Keeping  track  of  the  tangential  contact  displace¬ 
ment  history  emerges  as  a  key  element  of  the  model.  We 
show  that  identical  simulations  using  contact  models  that  in¬ 
clude  either  no  tangential  contact  displacement  history  or 
only  single -time -step  tangential  contact  displacement  history 
are  unable  to  accurately  model  the  direct  shear  test. 

1  The  Discrete  Element  Method 

Two  alternative  approaches  have  emerged  as  viable  so¬ 
lutions  for  large  frictional  contact  problems  in  granular  flow 
dynamics  and  quasi-static  geomechanics  applications.  The 
so-called  Complementarity  Method  (CM)  is  generally  fa¬ 
vored  within  the  multibody  dynamics  community,  see  for  in¬ 
stance  [1].  In  this  approach,  individual  particles  in  a  bulk 
granular  material  are  modeled  as  rigid  bodies,  and  non¬ 
penetration  conditions  are  written  as  complementarity  equa¬ 
tions  which,  in  conjunction  with  a  Coulomb  friction  law,  lead 
to  a  Differential  Variational  Inequality  (DVI)  form  of  the 
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Newton-Euler  equations  of  motion  [2] .  Not  limited  by  sta¬ 
bility  considerations,  DEM-CM  allows  for  much  larger  time 
integration  steps  than  the  alternative  Penalty  Method-based 
(PM)  solutions,  since  the  latter  involve  large  contact  stiff¬ 
nesses  that  impose  strict  stability  conditions  on  all  explicit 
time  integration  algorithms.  However,  DEM-CM  involves  a 
relatively  complex  and  computationally  costly  solution  se¬ 
quence  per  time  step,  since  it  leads  to  a  mathematical  pro¬ 
gram  with  complementarity  and  equality  constraints,  which 
must  be  relaxed  to  obtain  tractable  linear  complementarity  or 
cone  complementarity  problems  [3]. 

More  mature  and  widely  adopted  within  the  geomechan¬ 
ics  community  [4],  DEM-PM  can  be  viewed  either  as  a  reg¬ 
ularization  (or  smoothing)  approach,  which  relies  on  a  relax¬ 
ation  of  the  rigid-body  assumption,  or  as  a  deformable-body 
approach  localized  to  the  points  of  contact  between  individ¬ 
ual  particles  in  a  bulk  granular  material  [5,  6].  In  this  ap¬ 
proach,  normal  and  tangential  contact  forces  are  calculated 
using  various  laws  [7, 8],  which  are  based  on  the  local  body 
deformation  at  the  point  of  contact.  In  the  contact-normal 
direction,  this  local  body  deformation  is  defined  as  the  pen¬ 
etration  (overlap)  of  the  two  quasi-rigid  bodies.  In  the  tan¬ 
gential  direction,  the  deformation  is  defined  as  the  total  tan¬ 
gential  displacement  incurred  since  the  initiation  of  contact. 
Once  contact  forces  are  known,  the  time  evolution  of  each 
body  in  the  system  is  obtained  by  integrating  the  Newton- 
Euler  equations  of  motion.  Since  in  this  approach  the  con¬ 
tact  force-displacement  law  is  derived  from  the  elastic  prop¬ 
erties  of  the  materials  constituting  the  contacting  bodies;  i.e.. 
Young’s  modulus  and  Poisson’s  ratio,  the  DEM-PM  is  capa¬ 
ble  of  resolving  statically  indeterminate  loading  conditions 
that  can  exist  at  the  particle  level  in  random  granular  pack¬ 
ings  [9-11].  However,  due  to  large  contact  stiffnesses  and 
the  use  of  explicit  time  integration  [12],  the  DEM-PM  ap¬ 
proach  is  limited  to  very  small  time  integration  step- sizes  to 
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2  The  Penalty  Method  or  Soft-Body  Approach 

A  granular  or  particulate  medium  problem  is  modeled 
in  DEM  using  a  massive  collection  of  distinct  rigid  or  de¬ 
formable  elements  having  simple  shapes  that  in  many  cases 
are  spheres.  In  the  DEM-PM  or  soft-body  approach,  the  el¬ 
ements  are  “soft”  -  they  are  allowed  to  “overlap”  or  experi¬ 
ence  local  deformation  before  a  corrective  contact  force  is 
applied  at  the  point  of  contact.  Once  such  an  overlap 
is  detected,  by  any  one  of  a  number  of  contact  algorithms, 
contact  force  vectors  and  normal  and  tangential  to  the 
contact  plane  at  the  point  of  contact  are  calculated  using  var¬ 
ious  constitutive  laws  [7,  8]  based  on  the  local  body  defor¬ 
mation  at  the  point  of  contact.  In  the  contact-normal  direc¬ 
tion,  n,  this  local  body  deformation  is  defined  as  the  penetra¬ 
tion  (overlap)  of  the  two  quasi-rigid  bodies,  =  5„n.  In  the 
contact-tangential  direction,  the  deformation  is  defined  as  the 
total  tangential  displacement  incurred  since  the  initiation  of 
contact,  which  is  approximated  as  a  vector  in  the  contact 
plane. 

An  example  of  a  DEM-PM  contact  constitutive  law,  a 
slightly  modified  form  of  which  is  used  in  the  open-source 
codes  Chrono  [13]  and  LIGGGHTS  [14],  is  the  following 
viscoelastic  model  based  on  either  Hookean  or  Hertzian  con¬ 
tact  theory: 

^ n  —  f  {J^n^n  yn^^n) 

F(  =  f{R,8n)  {-k,u,  -y,m\t) , 

where  u  =  is  the  overlap  or  local  contact  displace¬ 

ment  of  two  interacting  bodies,  see  Pig.  1.  The  quantities 
fh  =  miJiij / {rui  +  nij)  and  R  =  RiRj / {Ri  +  Rj)  represent  the 
effective  mass  and  effective  radius  of  curvature,  respectively, 
for  contacting  bodies  with  masses  m/  and  mj  and  contact 
radii  of  curvature  Ri  and  Rj.  The  relative  velocity  at  the  con¬ 
tact  point,  V  =  +  v^,  and  its  normal  and  tangential  compo¬ 

nents  \n  and  v^  are  computed  as 

V  =  Vj  +  Q.j  X  Yj  —  V;  —  X  Yi 

v„  =  (n-v)n  (2) 

V/  =  V  -  v„ , 

where  V/  and  \j  are  the  velocity  vectors  of  the  centers  of 
mass  of  bodies  i  and  7,  Q^t  and  Qj  are  the  angular  velocity 
vectors  of  bodies  i  and  7,  and  and  Yj  are  the  position  vec¬ 
tors  from  the  centers  of  mass  of  bodies  i  and  7  to  the  point  of 
contact.  Por  Hookean  contact,  f{R,dn)  =  1  in  Eqn.  (1);  for 
Hertzian  contact,  one  can  let  /(^,5„)  =  [15, 16].  The 

normal  and  tangential  stiffness  and  damping  coefficients  kn, 
and  are  obtained,  through  various  constitutive  laws 
derived  from  contact  mechanics,  from  physically  measurable 
properties  for  the  materials  of  the  contacting  bodies,  such  as 
Young’s  modulus,  Poisson’s  ratio,  the  coefficient  of  restitu¬ 


tion,  etc.  Detailed  descriptions  of  the  contact  models  imple¬ 
mented  in  Chrono  and  LIGGGHTS,  as  well  as  alternative 
contact  models  are  provided  in  [17]. 

The  component  of  the  overlap  or  contact  displacement 
vector  u  in  the  contact-normal  direction,  =  5„n,  is  ob¬ 
tained  directly  from  the  contact  detection  algorithm,  which 
provides  the  magnitude  of  the  “inter-penetration”  5„.  It  fol¬ 
lows  that  Ufi  is  parallel  to  the  normal  component  of  the  rel¬ 
ative  velocity  vector  at  the  point  of  contact.  However,  it 
is  important  to  note  that,  in  general,  the  same  is  not  true  of 
the  tangential  component  of  the  overlap  vector,  or  tangential 
contact  displacement,  and  the  tangential  component  of  the 
relative  velocity  vector  .  Specifically,  they  must  lie  in  the 
contact  plane,  but  may  or  may  not  be  parallel  to  each  other. 
In  particular,  even  if  there  is  no  relative  tangential  velocity  at 
the  contact  point,  there  may  still  be  a  tangential  contact  force 
induced  by  deformation  in  the  tangential  plane,  and  this  force 
may  be  needed  to  support  static  friction. 

Herein,  the  tangential  contact  displacement  vector  is 
formulated  as 


where  t  is  the  current  time  and  is  the  time  at  the  initiation 
of  contact.  Por  the  true  tangential  contact  displacement  his¬ 
tory  model,  the  vector  must  be  stored  and  updated  at  each 
time  step  for  each  contact  point  on  a  given  pair  of  contacting 
bodies  from  the  time  that  contact  is  initiated  until  that  con¬ 
tact  is  broken.  The  tangential  (or  shear)  contact  displacement 
history  vector  is  then  given  at  time  step  i  by 

*  r  (4) 

where  Atf  is  the  integration  time  step  size,  tf  =  tt-i  -\-Ati,  and 
a  subscript  indicates  the  time  step  at  which  each  variable  is 
evaluated.  The  projection  of  u*  onto  the  contact  plane  is  nec¬ 
essary  to  ensure  that  is  in  the  contact  plane  at  each  time 
step.  Note  that  is  set  to  zero  at  the  initiation  of  contact, 
for  some  k. 

A  simpler  but  less  effective  tangential  contact  displace¬ 
ment  model  suggested  in  the  literature  is  a  single  time  step 
approximation  of  Eqn.  (4),  given  for  any  time  step  by 

Ut  =  \tAt.  (5) 

This  model,  which  we  will  call  pseudo-history,  essentially 
assumes  that  contact  never  persists  for  more  than  a  single 
time  step,  and  it  is  unable  to  support  a  static  friction  force  in 
the  absence  of  relative  tangential  velocity. 

To  enforce  the  Coulomb  friction  law,  if  |F,|  >/i|F„|  at 
any  given  time  step,  then  before  the  contributions  of  the  con¬ 
tact  forces  are  added  to  the  resultant  force  and  torque  on  the 
body,  the  (stored)  value  of  |u^|  is  scaled  so  that  |F^|  =7/|F„|, 
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Fig.  1.  DEM-PM  contact  model  described  in  this  section,  with  nor¬ 
mal  overlap  distance  5„,  contact-normal  unit  vector  n,  and  tan¬ 
gential  displacement  vector  in  the  plane  of  contact  (top),  and 
with  a  Hookean-linear  contact  force-displacement  law  with  constant 
Coulomb  sliding  friction  (bottom). 


where  ju  is  the  Coulomb  (static  and  sliding)  friction  coeffi¬ 
cient.  For  example,  if  f{x)  =  1  in  Eqn.  (1),  then 


- LIGGGHTS:  True  History 


kt\\lt\>  iu\Fn\  ^  .  (6) 

Figure  1  illustrates  the  DEM-PM  contact  model  de¬ 
scribed  in  this  section  with  a  Hookean-linear  contact  force- 
displacement  law  with  constant  Coulomb  sliding  friction. 
Once  the  contact  forces  and  are  computed  for  each 
contact  and  their  contributions  are  summed  to  obtain  a  resul¬ 
tant  force  and  torque  on  each  body  in  the  system,  the  time 
evolution  of  each  body  in  the  system  is  obtained  by  inte¬ 
grating  the  Newton-Euler  equations  of  motion,  subject  to  the 
Courant-Eriedrichs-Lewy  (CEL)  stability  condition,  which 
requires  [18]  that  At  <  A^crit  ^  \/ ^min/ ^max- 

3  The  Importance  of  Multi-Step  Tangential  Contact 

Displacement  History 

To  demonstrate  the  importance  of  using  tangential  dis¬ 
placement  history  in  the  DEM-PM  contact  model,  we 
first  perform  direct  shear  simulations  of  small  randomly 
packed  specimens  of  1,800  and  5,000  identical  spheres  in 
Chrono  [13]  and  LIGGGHTS  [14].  The  inside  dimensions 
of  the  shear  box  are  6  cm  in  length  by  6  cm  in  width,  and 
the  height  of  the  granular  material  specimen  is  also  approxi¬ 
mately  6  cm.  The  spheres  have  a  uniform  diameter  of  5  mm. 
The  random  packing  of  1 , 800  spheres  was  initially  obtained 
by  a  “rainfall”  method,  after  which  the  spheres  were  com¬ 
pacted  with  friction  temporarily  turned  off  to  obtain  a  dense 
packing.  The  resulting  void  ratio  was  approximately  e  =  0.4, 
which  corresponds  to  a  dense  packing  [19,20].  Eor  this  com¬ 
parison,  the  material  properties  for  spheres  were  taken  to  be 
those  corresponding  to  quartz  -  the  density  is  2,500  kg/no?, 


Fig.  2.  Direct  shear  simulation  setup  (top)  and  shear  versus 
displacement  results  (bottom)  obtained  by  Chrono  [13]  and 
LIGGGHTS  [14]  for  1,800  randomly  packed  uniform  spheres 
using  the  true  tangential  contact  displacement  history  model  of 
Eqn.  (4),  the  pseudo-history  model  of  Eqn.  (5),  and  no  tangential 
contact  history. 


the  inter-particle  friction  coefficient  is  jli  =  0.5,  Poisson’s  ra¬ 
tio  is  V  =  0.3,  and  the  elastic  modulus  is  E  =  8(10^^)  Pa. 
However,  in  order  to  ensure  a  stable  simulation  with  a  rea¬ 
sonable  time  integration  step-size  of  At  =  10“^  s,  the  elas¬ 
tic  modulus  was  reduced  by  four  orders  of  magnitude  to 
E  =  8(10^)  Pa.  The  shear  speed  was  1  mm/s.  The  simu¬ 
lation  geometry  in  its  final  position  is  shown  in  Pig.  2  (top). 
Pigure  2  (bottom)  shows  the  shear-displacement  curves  ob¬ 
tained  by  Chrono  and  LIGGGHTS  with  the  same  normal 
and  tangential  contact  force-displacement  models.  The  la¬ 
bels  “True  History”  and  “No  History”  refer  to  whether  or  not 
tangential  contact  history  is  stored  and  used  in  the  friction 
model.  Included  in  this  comparison  is  the  “Pseudo-History” 
scenario,  in  which  the  tangential  contact  displacement  vec¬ 
tor  is  approximated  by  the  product  of  the  relative  tangential 
velocity  vector  at  the  contact  point  and  the  time  step- size  at 
any  given  time.  This  pseudo-history  approach  is  attractive, 
since  unlike  the  “True  History”  alternative,  it  avoids  the  stor¬ 
age  of  a  tangential  contact  history  vector  over  multiple  time 
steps  for  each  contact  point.  However,  Pig.  2  shows  that  the 
pseudo-history  approximation  is  no  better  than  ignoring  tan¬ 
gential  displacement  history  altogether  for  the  quasi- static 
direct  shear  test.  This  is  explained  by  the  observation  that 
under  quasi-static  (or  static)  deformation  conditions,  the  de- 
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pendence  of  the  pseudo-history  approximation  on  the  rela¬ 
tive  inter-particle  tangential  velocity,  which  is  zero,  effec¬ 
tively  eliminates  the  inter-particle  tangential  contact  force 
and  so  renders  the  inter-particle  friction  coefficient  ju  effec¬ 
tively  zero. 

Also  noteworthy  in  Fig.  2  is  the  fact  that  the  inter¬ 
particle  friction  coefficient  /u  for  the  spheres,  which  could 
also  be  described  as  a  micro-scale  “inter-particle  friction 
angle”  (|)^  =  tan“^^,  is  nowhere  close  to  having  the  same 
value  as  the  macro-scale  “material  friction  coefficient”  ^macro 
for  the  bulk  granular  material.  The  latter,  more  commonly 
described  as  a  bulk  granular  material  friction  angle  (|)  = 
tan“^/iniacro.  is  the  material  parameter  that  defines  the  yield 
surface  for  the  bulk  granular  material  according  to  the  Mohr- 
Coulomb  yield  criterion.  The  material  friction  angle  (|)  is 
also  known  as  the  angle  of  repose  for  the  bulk  granular  ma¬ 
terial.  Nor  should  it  be  surprising  that  (|)  /  (|)^,  since,  as 
noted  in  [21],  even  if  the  inter-particle  friction  coefficient 
jLi,  and  hence  the  micro-scale  friction  angle  (|)^,  is  zero,  the 
bulk  granular  material  friction  angle  (|)  will  in  general  not  be 
zero.  Rather,  if  ^  =  0,  then  (|)  =  tj/,  where  \\f  is  the  dilation  an¬ 
gle  of  the  granular  material.  Note  that  typically,  ij/  ^  15° 
for  densely  packed  well-graded  sands  [22].  In  particular, 
we  note  from  Fig.  2  that,  when  the  tangential  contact  dis¬ 
placement  history  model  is  used,  while  ju  =  0.5  and  hence 
(|)^  ^26.6°  for  the  spheres,  the  peak  ratio  of  shear  stress  to 
normal  stress  for  the  bulk  granular  material  is  ^macro  ~  2,  and 
hence  (^p  ^  63°;  and  the  residual  ratio  of  shear  stress  to  nor¬ 
mal  stress  for  the  bulk  granular  material  is  ^macro  ~  and 
hence  (|)r  ~  45°.  On  the  other  hand,  when  the  tangential  con¬ 
tact  displacement  history  model  is  not  used,  ^macro  ~  0.25 
throughout  the  simulation,  and  hence  (|)^  =  (|)^  =  (|)  ^  14°. 
Note  that  all  of  these  results  are  obtained  in  the  absence  of 
any  rolling  or  spinning  friction. 

To  emphasize  the  importance  of  using  multi-step  tan¬ 
gential  contact  displacement  history,  it  should  be  pointed  out 
that  other  factors  involved  in  the  model,  such  as  the  values 
of  kn,  kf,  Jn,  and  y^,  turned  out  to  play  a  secondary  role  in 
the  outcome  of  the  simulation.  In  fact,  a  significant  degree 
of  variation  exists  in  the  literature  for  the  exact  values  of  the 
contact  stiffness  coefficients  kn  and  kt  [17].  The  same  is  true 
for  the  mass  proportional  damping  coefficients  y„  and  Jt  •  The 
latter  are  frequently  simply  chosen  sufficiently  large  to  elim¬ 
inate  numerical  noise  in  the  DEM-PM  simulations.  For  ex¬ 
ample,  the  results  of  DEM-PM  simulations  of  direct  (ring) 
shear  tests  with  periodic  boundary  conditions  on  ASTM  C 
778-06  standard  graded  (quartz)  sand  with  a  log-normal  par¬ 
ticle  size  distribution,  mean  diameter  £>50  =  0.35  mm,  and 
coefficient  of  uniformity  Cu  =  1.7  were  considered  in  [23]. 
In  these  simulations,  which  employed  the  multi-step  tangen¬ 
tial  contact  displacement  history  model  described  herein,  the 
damping  coefficients  in  Eqn.  (1)  were  taken  to  be  y„  =  40  s“^ 
and  jt  =20  s“\  and  the  contact  stiffnesses  kn  and  kt  were 
taken  to  be  constant,  with  kn  =  10^^  N/m  and  kt  =  8(10^^) 
N/m.  Despite  these  simplifications,  and  the  fact  that  the  sim¬ 
ulations  performed  included  no  rolling  friction  and  the  sand 
particles  were  modeled  as  spheres  of  different  sizes,  the  cor¬ 
rect  macro-scale  residual  bulk  granular  material  friction  an¬ 


gle  of  (|)r  =  30°  [24]  was  reproduced  exactly.  The  only  other 
material  parameter  that  needed  to  be  specified,  in  addition 
to  the  particle  size  distribution,  was  the  inter-particle  fric¬ 
tion  coefficient  ^  =  0.5,  which  is  considered  by  Mitchell  and 
Soga  [20]  to  be  “reasonable  for  quartz,  both  wet  and  dry.” 
Note  that  the  values  of  the  peak  and  residual  friction  angles 
are  strongly  dependent  on  the  particle  size  distribution  [25], 
which  is  why  the  residual  friction  angle  for  uniform  quartz 
spheres  cannot  be  expected  to  be  the  same  as  that  of  quartz 
spheres  (or  well-rounded  quartz  sand)  with  a  log-normal  par¬ 
ticle  size  distribution. 


4  Validation  Against  Direct  Shear  Experiments  With 

Uniform  Glass  Beads 

Whereas  the  previous  section  demonstrated  the  differ¬ 
ence  in  results  between  the  “True  History”  and  “No  His¬ 
tory”  scenarios,  herein  we  compare  the  “True  History” 
shear-displacement  curves  against  experimental  data  re¬ 
ported  in  [26].  Specifically,  to  verify  that  the  Chrono 
DEM-PM  contact  model  with  true  tangential  displacement 
history  currently  does  indeed  accurately  model  the  micro¬ 
scale  physics  and  emergent  macro-scale  properties  of  a  sim¬ 
ple  granular  material.  Fig.  3  shows  shear  versus  displace¬ 
ment  curves  obtained  from  both  experimental  [26]  (top)  and 
Chro no-simulated  (center  and  bottom)  direct  shear  tests, 
performed  under  constant  normal  stresses  of  3.1,  6.4,  12.5, 
and  24.2  kPa,  on  5,000  uniform  glass  beads.  The  simulation 
geometry  in  its  final  position  is  similar  to  that  shown  in  Fig.  2 
(top),  except  that  the  inside  dimensions  of  the  shear  box  are 
now  12  cm  in  length  by  12  cm  in  width.  The  height  of  the 
granular  material  specimen  in  the  box  is  still  approximately 
6  cm.  In  both  the  experimental  and  simulated  direct  shear 
tests,  the  glass  spheres  have  a  uniform  diameter  of  6  mm, 
and  the  random  packing  of  5,000  spheres  was  initially  ob¬ 
tained  by  a  “rainfall”  method,  after  which  the  spheres  were 
compacted  by  the  confining  normal  stress  without  adjusting 
the  inter-particle  friction  coefficient.  The  DEM-PM  simu¬ 
lations  were  performed  in  Chrono  using  a  Hertzian  normal 
contact  force  model  and  true  tangential  contact  displacement 
history  with  Coulomb  friction.  The  material  properties  of  the 
spheres  in  the  simulations  were  taken  to  be  those  correspond¬ 
ing  to  glass  [26],  for  which  the  density  is  2,550  kg/m^,  the 
inter-particle  friction  coefficient  is  /i  =  0.18,  Poisson’s  ratio 
is  V  =  0.22,  and  the  elastic  modulus  is  E  =  4(10^^)  Pa,  except 
that  the  elastic  modulus  was  again  reduced  by  several  orders 
of  magnitude,  to  E  =  4(10^)  Pa  (center)  and  E  =  4(10^)  Pa 
(bottom)  for  comparison,  to  ensure  a  stable  simulation  with 
a  reasonable  time  integration  step-size  of  Ar  =  10“^  s.  The 
shear  speed  was  1  mm/s. 

Eigure  3  (center)  shows  that  the  DEM-PM  direct  shear 
simulations  performed  in  Chrono  on  5,000  glass  spheres 
with  E  =  4(10^)  Pa  matches  reasonably  well  the  physical 
experiments  for  all  but  the  highest  normal  stress  of  24.2  kPa. 
This  observed  error  in  the  simulation  results,  which  increases 
with  increasing  normal  stress,  is  consistent  with  the  fact  that 
the  contact  stiffness  for  the  spheres  in  these  DEM-PM  simu¬ 
lations  is  four  orders  of  magnitude  smaller  than  the  stiffness 
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correct  contact  stiffness,  results  in  a  peak  and  residual  shear 
stress  that  is  much  closer  to  the  experimentally  observed  val¬ 
ues  for  all  four  of  the  constant  normal  stresses  tested.  This 
is  a  significant  observation,  since  it  has  often  been  argued  in 
the  DEM-PM  literature  that  decreasing  the  value  of  the  elas¬ 
tic  modulus  to  allow  a  larger  stable  time  step-size  should  only 
affect  the  elastic  portion  of  the  shear  displacement  curve  for 
the  bulk  granular  material.  A  comparison  of  Figs.  3  (center) 
and  3  (bottom),  however,  while  confirming  this  difference  in 
the  elastic  portion  of  the  shear-displacement  curve,  also  re¬ 
veals  a  significant  difference  in  the  plastic  or  post-yield  por¬ 
tion  of  the  shear-displacement  curve  for  the  direct  shear  test, 
in  particular  the  peak  and  residual  shear  stresses,  and  the  cor¬ 
responding  peak  and  residual  friction  angles,  for  all  four  of 
the  constant  normal  stresses  tested. 


5  Conclusions 

In  relation  to  using  computer  simulation  to  capture  the 
dynamics  of  granular  material,  this  technical  note  makes  the 
following  two  points.  First  and  foremost,  contrary  to  com¬ 
mon  perception,  eliminating  the  tangential  contact  history 
in  DEM-PM  yields  wrong  results  in  a  shear  test  that,  while 
basic  and  deceptively  simple,  remains  difficult  to  simulate. 
Moreover,  a  quasi-history  approach  that  only  relies  on  the 
tangential  deformation  at  the  current  time  step  produces  in¬ 
accurate  results.  Second,  a  comparison  against  experimental 
data  suggests  that  the  simulation  results  are  only  moderately 
impacted  by  the  values  selected  for  the  DEM-PM  model  pa¬ 
rameters,  of  which  the  normal  stiffness  kn  turns  out  to  quan¬ 
titatively  infiuence  the  most  the  outcome  of  the  numerical 
experiments.  Specifically,  over  a  broad  spectrum  of  values 
for  y„,  and  y^,  the  simulation  results  are  qualitatively  ac¬ 
ceptable  for  artificially  low  values  of  a  compromise  made 
in  order  to  allow  stable  numerical  integration  at  larger  sim¬ 
ulation  time  steps.  However,  more  accurate  results  call  for 
higher  values  of  kn  that  come  close  to  the  theory  predicted 
values  for  this  parameter. 


Fig.  3.  Direct  shear  test  results  for  5,000  randomly  packed  uniform 
glass  beads  obtained  by  experiment  [26]  (top)  and  DEM-PM  simu¬ 
lations  using  Chrono  (center  and  bottom),  under  constant  normal 
stresses  of  3.1,  6.4,  12.5,  and  24.2  kPa.  For  the  DEM-PM  simula¬ 
tions,  elastic  moduli  of  £”  =  4(10^)  Pa  (center)  and  E  =  4(10^)  Pa 
(bottom)  are  used. 
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of  true  glass  beads.  To  explore  the  effect  that  the  value  of  the 
elastic  modulus  has  on  the  DEM-PM  direct  shear  results,  we 
have  also  performed  DEM-PM  simulations  using  an  elastic 
modulus  of  £■  =  4(  10^)  Pa  for  the  spheres,  which  is  still  three 
orders  of  magnitude  smaller  than  the  true  elastic  modulus 
of  glass  beads.  Figure  3  (bottom)  shows  that  increasing  the 
value  of  the  elastic  modulus  of  the  spheres  in  the  direct  shear 
simulations  by  an  order  of  magnitude  to  £"  =  4(10^)  Pa;  i.e., 
using  a  contact  stiffness  for  the  spheres  that  is  three  rather 
than  four  orders  of  magnitude  smaller  than  the  physically 
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